HUTP-05/A0038 
MIT-CTP 3670 



X 

5-S 



,a,c 



Limits on non-Gaussianities from WMAP data 

\^ . Paolo Creminelli a , Alberto Nicolis a , Leonardo Senatore b,d 

q ; Max Tegmark d and Matias Zaldarriaga 8 
(N 

Q_i' a Jefferson Physical Laboratory, 

■ Harvard University, Cambridge, MA 02138, USA 

b Center for Theoretical Physics, 
Massachusetts Institute of Technology, Cambridge, MA 02139, USA 

> 

Q\ ' c Center for Astrophysics, 

(N ' 

o i 

d Department of Physics, 

■ Massachusetts Institute of Technology, Cambridge, MA 02139, USA 

O 

Oh; 

O ' Abstract 



We develop a method to constrain the level of non-Gaussianity of density perturbations when the 
3-point function is of the "equilateral" type. Departures from Gaussianity of this form are produced 



by single field models such as ghost or DBI inflation and in general by the presence of higher order 
derivative operators in the effective Lagrangian of the inflaton. We show that the induced shape 
of the 3-point function can be very well approximated by a factorizable form, making the analysis 
practical. We also show that, unless one has a full sky map with uniform noise, in order to saturate 
the Cramer-Rao bound for the error on the amplitude of the 3-point function, the estimator must 
contain a piece that is linear in the data. We apply our technique to the WMAP data obtaining 
a constraint on the amplitude f^ lL of "equilateral" non-Gaussianity: —366 < /ni!" 1 ^ a ^ 
95% C.L. We also apply our technique to constrain the so-called "local" shape, which is predicted 
for example by the curvaton and variable decay width models. We show that the inclusion of the 
linear piece in the estimator improves the constraint over those obtained by the WMAP team, to 
-27 < /^; al < 121 at 95% C.L. 



1 Introduction 

Despite major improvements in the quality of cosmological data over the last few years, there are very 
few observables that can characterize the early phases of the Universe, when density perturbations 
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were generated. Deviations from a purely Gaussian statistics of density perturbations would be a 
very important constraint on models of early cosmology. In single field slow-roll inflation, the level 
of non-Gaussianity is sharply predicted ^ Hj to be very small, less than This is quite far 

from the present experimental sensitivity. On the other hand, many models have recently been 
proposed with a much higher level of non-Gaussianity, within reach of present or forthcoming data. 
For nearly Gaussian fluctuations, the quantity most sensitive to departures from perfect Gaussianity 
is the 3-point correlation function. In general, each model will give a different correlation between 
the Newtonian potential modes 1 : 

($(ki)$(k 2 )$(k 3 )> = (2ir) 3 5 3 (k 1 + k 2 + k 3 )F(k 1 ,k 2 ,k 3 ) . (1) 

The function F describes the correlation as a function of the triangle shape in momentum space. 

The predictions for the function F in different models divide quite sharply into two qualitatively 
different classes as a consequence of qualitatively different ways of producing correlations among 
modes [4_. The first possibility is that the source of density perturbations is not the inflaton but a 
second light scalar field a. In this case non-Gaussianities are generated by the non-linear relation 
between the fluctuation 5a of this field and the final perturbation $ we observe. This non-linearity 
is local as it acts when the modes are much outside the horizon; schematically we have 3?(x) = 
A5o~(x) + -B(<5ct 2 (x) — (5a 2 )) + . . . . Even starting from a purely Gaussian 8a, the quadratic piece 
introduces a 3-point function for $ of the form 

F(h, k 2 , h) = ■ 2A| • (J, + J, + JL) , (2) 

where A$ is the power spectrum normalization, (<E , (ki)$(k2)) = (27r) 3 <5 3 (ki + k 2 )A$ • /cj~ 3 , which 
for the moment has been taken as exactly scale invariant, and where /^l^ 1 is proportional to B. 
Examples of this mechanism are the curvaton scenario and the variable decay width model 0, 
which naturally give rise to /nl^ greater than 10 and 5, respectively. 

The second class of models are single field models with a non-minimal Lagrangian, where the 
correlation among modes is created by higher derivative operators E3 E3 ■ I n this case, the 
correlation is strong among modes with comparable wavelength and it decays when we take one of fc's 
to zero keeping the other two fixed. Although different models of this class give a different function 
F, all these functions are qualitatively very similar. We will call this kind of functions equilateral: 
as we will see, the signal is maximal for equilateral configurations in Fourier space, whereas for the 
local form (J2J) the most relevant configurations are the squeezed triangles with one side much smaller 
than the others. 

The strongest constraint on the level of non-Gaussianity comes from the WMAP experiment. 
The collaboration analyzed their data searching for non-Gaussianity of the local form (j2j), finding 
the data to be consistent with purely Gaussian statistics and placing limits on the parameter /^l^ 

m 

-58 < /^ al < 134 at 95% C.L. (3) 

^^Even with perfectly Gaussian primordial fluctuations, the observables, e.g. the temperature anisotropy, 
will not be perfectly Gaussian as a consequence of the non-linear relation between primordial perturbations 
and what we will eventually observe. These effects are usually of order 10 -5 (see for example and thus 
beyond (but not much) present sensitivity. In the following we will disregard these contributions. 
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The main purpose of this paper is to perform a similar analysis searching for non-Gaussianities 
of the equilateral form. We can extend the definition of /nl^ in eq. P|) to a generic function F by 
setting the overall normalization on equilateral configurations: 

6A 2 

F(k,k,k) = f NL - 1 f. (4) 

In this way, two different models with the same /nl will give the same 3-point function for equilateral 
configurations. For equilateral models in particular, the overall amplitude will be characterized by 
/jjT 111 '- Indirect constraints on /nl" 1 ' have been obtained starting from the limit of eq. © in [I], 
resulting in | /nl 111 I ~ 800; we will see that a dedicated analysis is, as expected, much more sensitive. 
Many of the equilateral models naturally predict sizeable values of the parameter f^jf*"'. ghost 
inflation and DBI inflation [315] tend to have /ni!" 1 ~ 100, tilted ghost inflation f^ lL > 1 [TO] , 
while slow roll inflation with higher derivative coupling typically give /nl^ ^ 1 Cl- 
in Section |21 we discuss the main idea of the analysis. We will see that an optimal analysis is 
numerically very challenging for a generic form of the function F, but simplifies dramatically if the 
function F is factorizable in a sense that will be defined below. As all the equilateral forms predicted 
in different models are (qualitatively and quantitatively) very similar, our approach will be to choose 
a factorizable function that well approximates this class. The analysis is further complicated by the 
breaking of rotational invariance: a portion of the sky is fully masked by the presence of the galaxy 
and moreover each point is observed a different number of times. In Section|3]we look for an optimal 
estimator for /nl- We discover that, unlike the rotationally invariant case, the minimum variance 
estimator contains not only terms which are cubic in the temperature fluctuations, but also linear 
pieces. These techniques are used to analyze the WMAP data in Section 0J The result is that the 
data are compatible with Gaussian statistics. The limit for the equilateral models is 



-366 < / N q L uil - < 238 at 95% C.L. (5) 



We also obtain a limit on /nl 3,1 



-27 < /^ al < 121 at 95% C.L. (6) 

which is a slight improvement with respect to (JSJ) as a consequence of the above-mentioned additional 
linear piece in the estimator. 



2 A factorizable equilateral shape 

It is in principle straightforward to generalize the analysis for non-Gaussianities from the local 
model to another one. We start assuming that the experimental noise on the temperature maps is 
isotropic and that the entire sky is observed (no Galactic and bright source mask); we will relax 
these assumptions in the next Section. In this case it can be proved that an optimal estimator £ for 
/nl exists |12l I13j , that is an estimator which saturates the Cramer- Rao inequality and thus gives 
the strongest possible constraint on the amplitude of the non-Gaussian signal. The estimator £ is a 
sum of terms cubic in the temperature fluctuations, each term weighted by its signal to noise ratio: 



-■y 



{a, 



km, 11 12 ' 3 



(7) 
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Given the assumptions, the power spectrum Q (which is the sum of the CMB signal and noise) is 
diagonal in Fourier space; TV" is a normalization factor which makes the estimator unbiased. From 
rotational invariance we can simplify the expression introducing the Wigner 3j symbols to 

P-l VI ' ! ' 2 ^3 | B hhh (9 \ 
N ^ \ mi m 2 m 3 / C h Ci 2 C h 

where Bi { i 2 i 3 is the angle-averaged bispectrum which contains all the information about the model 
of non-Gaussianity we are considering. If B^i 2 i 3 is calculated for /nl = 1, then the normalization 
factor N is given by 

We now have to relate the angle-averaged bispectrum to the underlying correlation among 3d modes 
F(ki,k 2 ,kz). After some manipulations following |14j . the estimator takes the form 

oo 

£ = Jf ^J d 2 nY hmi (n)Y hm2 (h)Y hm3 (h) J r 2 dr j h (hr)ji 2 (k 2 r)j h (k 3 r) C^C^C^ 



kn~ii 







2kjdh 2kldk2 2kldh F{ kl , k 2 , k 3 )Al(h)Al{k 2 )Al(k 3 ) a hmi a l2m2 a l3m3 , (10) 



7T 7T 7T 



where Aj (k) is the CMB transfer function which relates the a\ m to the Newtonian potential Q(k): 

^Af(kMk)YL(k) . (11) 

Equation (|1U() is valid for any shape of the 3-point function F. Unfortunately this expression is 
computationally very challenging. The sums over m can be taken inside the integrals and factorized, 
but we are still left with a triple sum over / of an integral over the sphere. The calculation time 
grows as N^ 2 els , where the number of pixels -/V p i xe i s is of order 3 x 10 6 for WMAP. This approach is 
therefore numerically too demanding. 

As noted in Jl], a crucial simplification is possible if the function F is factorizable as a product 
of functions of k\, k 2 and k% or can be written as a sum of a small number of terms with this property. 
In this case the second line of (|1U|) becomes factored as the product of functions of each / separately, 
so that now also the sum over / can be done before integrating over the sphere. For example if we 
assume that F{ki,k 2 ,k%) = f\{ki)f 2 {k 2 )f^{k^), the estimator simplifies to 



£ = -.j d 2 n J r 2 dr ]}£ J ^j^f^Alik^a^Y^n) . (12) 

i=l Urn. 



3 /2 

The calculation is obviously much faster now: it now scales like -/V" } xe j s and it is dominated by going 
back and forth between real and spherical harmonics space. From expression we see that this 
simplification is possible for the local shape, and it was indeed used for the analysis of the WMAP 
data in [Tlllin]- 

Unfortunately, none of the "equilateral models" discussed in the Introduction predicts a function 
F which is factorizable (see some explicit expressions in 4 ), so that it is not easy to perform an 
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optimal analysis for a particular given model. However, all these models give 3-point functions which 
are quite similar, so that it is a very good approximation to take a factorizable shape function F 
which is close to the class of functions we are interested in and perform the analysis for this shape. 
In the limit k\ — ► with k 2 and k% fixed, all the equilateral functions diverge as k^ 1 @] (while the 
local form eq. @ goes as &f 3 ). An example of a function which has this behavior, is symmetric in 
fei, k 2 and £3, and is a sum of factorizable functions (and is homogeneous of order k~ 6 , see below) 
is 2 " 

F(h, k 2 , k 3 ) = • 6A| • ("pL - pp ~ pi3 - + J-L-p + ( 5 Perm.) j , (14) 

\ ft-^ r\j<2 n, 1 ftg ft 2 'vg rvj ft 9 ftg ft- 1 ft 2 ' b 3 / 

where the permutations act only on the last term in parentheses. The mild divergence for k\ — > is 
achieved through a cancellation among the various terms. 

In figure ^ we compare this function with the local shape. The dependence of both functions 
under a common rescaling of all k's is fixed to be cx A; -6 by scale invariance, so that we can factor out 
fcj -6 for example. Everything will now depend only on the ratios k 2 /k\ and k^/k%, which fix the shape 
of the triangle in momentum space. For each shape we plot F(l, k 2 /k±, k 3 /ki){k 2 /k\) 2 {k 3 /ki) 2 ; this 
is the relevant quantity if we are interested in the relative importance of different triangular shapes. 
The square of this function in fact gives the signal to noise contribution of a particular shape in 
momentum space 4,. We see that for the function (|14|). the signal to noise is concentrated on 
equilateral configurations, while squeezed triangles with one side much smaller than the others are 
the most relevant for the local shape. 

In figure |2] we study the equilateral function predicted both in the presence of higher-derivative 
terms [7_ and in DBI inflation 9 . In the second part of the figure we show the difference between 
this function and the factorizable one used in our analysis. We see that the relative difference is 
quite small. The same remains true for other equilateral shapes (see jl] for the analogous plots for 
other models). 

In 4 , a "cosine" between different shapes was defined which quantifies how different is the signal 
given by two distributions. The cosine is calculated from the scalar product of the functions in fig-H 
We can think about this cosine as a sort of correlation coefficient: if the cosine is close to 1 the two 
shapes are very difficult to distinguish experimentally and an optimal analysis for one of them is 
very good also for the other. On the other hand, a small cosine means that, once non-Gaussianities 
are detected, there is a good chance to distinguish the two functions and that an optimal analysis for 
one shape is not very good for the other. The cosine between our template shape and the functions 
predicted by equilateral models is very close to one (0.98 with the ghost inflation [H] prediction and 

2 Eq. (|14|) can be derived as follows. In order to make the divergence of F mild in the squeezed limit we 
can use at the numerator a quantity which goes to zero in the same limit. The area of the triangle does the 
job, going like k\ for k\ — > 0. The area can be expressed purely in terms of the sides through Heron's formula 
|16| . A = y/ s(s — k\){s — k2){s — k-i), where s = \{k\ + k 2 + k^) is the semiperimeter. The first s in the 
square root is irrelevant for our purposes, since it goes to a constant in the squeezed limit; we will therefore 
omit it. Also we want a sum of factorizable functions, so we get rid of the square root by considering A 2 . In 
conclusion, a function with all the features stated above is 

j?(u u u\ (s - ki)(s - k 2 )(s - fc 3 ) , . 

F(ki, k 2 ,k 3 ) cx , (13) 

"a "-2 ^3 

which, once expanded, reduces exactly to eq. (|14f> . 
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Figure 1: Plot of the function F(l, fe/fci, &3 Ai)(^2 Ai) 2 (&3Ai) 2 for the equilateral shape used in 
the analysis (top) and for the local shape (bottom). The functions are both normalized to unity for 
equilateral configurations |^ = |^ = 1. Since F(k\, &2, ks) is symmetric in its three arguments, it 
is sufficient to specify it for k\ > k<i > &3, so |^ < |^ < 1 above. Moreover, the triangle inequality 
says that no side can be longer than the sum of the other two, so we only plot F in the triangular 
region 1 — ^ < y- < < 1 above, setting it to zero elsewhere. 
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Figure 2: Top. Plot of the function F(l, fa/fci, k^/ki)(k2/ki) 2 {k^/ki) 2 predicted by the higher- 
derivative |.7~ and the DBI models Bottom. Difference between the above plot and the analogous 
one (top of fig. for the factorizable equilateral shape used in the analysis. 
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0.99 for higher derivative/DBI models 005]). This means that the error introduced in the analysis 
by the use of the factorizable shape instead of the correct prediction for a given equilateral model 
is at the percent level. On the other hand, as evident from figure Q ° ur template shape is quite 
different from the local model — the cosine is merely 0.41. 

All these numbers are obtained in 3 dimensions assuming that we can directly measure the 
fluctuations with a 3d experiment like a galaxy survey. The CMB anisotropics are a complicated 
projection from the 3d modes. This makes it more difficult to distinguish different shapes, although 
the picture remains qualitatively the same. 

Let us proceed with the study of our template shape. To further simplify the estimator in 
eq. pup, we define the functions 

ai(r) = - / dkk 2 Af(k)j l (kr) (15) 

7T JO 

A(r) = -/ dkk- 1 Af(k)j l (kr)A 9 (16) 



7T 



2 



o 



7T 



Ti(r) = -/ dkkAfik^krJA 1 ^ (17) 



o 



9 r+oo 

<y,(r) = - / dkAj(k)j l (kr)A 2 ^ . (18) 

We use this strange ordering of the functions to keep the notation compatible with [111 I15j . where 
the functions a and (3 were introduced for the analysis of the local shape. To evaluate £, we start 
from the spherical harmonic coefficients of the map ai m and calculate the four maps 

Mr, n) = Y, ^P-Ylm(h)a lm , B(r, n) = £ ^Y lm (n)a lm , (19) 

Ira Im 

C(r, n) = J2 ^Y lm {h)a lm , D(r, n) = ^ 6 -^-Y lm (h)a lm . (20) 

Ira Ira 

Now the estimator £ is given by 



£ = — — I r 2 dr I d 2 h 



2 

. ... .. .. A(r,n)B(r,h) 2 + -D(r,h) 3 - 2B(r,h)C(r,h)D(r,h) 

iV J J 3 

and the normalization N can be calculated from Q using the explicit form 



(21) 



Bhhh _ v - (Ml + i)(^ 1)W , + i) ( , , k j x (22) 

/>oo 

6 / r 2 dr [-a h (r)/3 h (r)(3 h (r) + (2 perm.) - 2<5j 1 (r)<5/ 2 (r)<5; 3 (r) + Ph{r)~ii 2 {r)5 h (r) + (5 perm.)] . 

JO 

The functions a, /3, 7 and 5 used in the analysis can be obtained numerically starting from the 
transfer functions Af(k), which can be computed given a particular cosmological model with publicly 
available software as CMBFAST 17 j. Some plots of these functions are given in fig. |5J where we 
choose values of r close to tq — tr (conformal time difference between recombination and the present), 
as these give the largest contribution to the estimator. The oscillatory behavior induced by the 
transfer function A[(k) is evident, with a peak at / ~ 200. The factors of 1(1 + 1) and (21 + 1) in 





8 / - 3 f r -0.4r R 

Figure 3: The functions a;(r) (in units of Mpc~ 3 ), (3i(r) (dimensionless), ji(r) (in units of Mpc~ 2 ), 
and 5i(r) (in units of Mpc -1 ) are shown for various radii r as functions of the multipole number I. 
The cosmological parameters are the same ones used in the analysis: f^/i 2 = 0.024, £l m h 2 = 0.14, 
h = 0.72, r = 0.17. With these parameters, the present conformal time To is 13.24 Gpc and the 
recombination time tr is 0.27 Gpc. 



the f3 and 5 functions, respectively, can be understood from the behavior at low I's, in the Sachs- 
Wolfe regime. Here the transfer function can be approximated by A-f(k) = —ji (fc(ro — tr)) /3, 
and we then see that (3i(tq — tr) oc f^°dk k~ l jf{k{jo — tr)) oc 1/ (/(/ + 1)) and <5/(to — tr) oc 
f^°dk jf(k(ro — tr)) oc 1/(2/ + 1). From these expressions we also see that the function (3 (the 
only one which is dimensionless) is very similar to the C; cmt "s as, in the Sachs- Wolfe regime, Cf mh ~ 

^J^dkk-l A,J?(fc(7b-7*)). 

3 Optimal estimator for / NL 

There are two quite general experimental complications that make the estimator £ defined in the last 
Section non-optimal. First, foreground emission from the Galactic plane and from isolated bright 
sources contaminates a substantial fraction of sky, which must be masked out before the analysis |18j . 
This projection which can be accomplished by giving very large noise to regions that are affected 
by foregrounds, besides reducing the available amount of data, has the important effect of breaking 
rotational invariance, so that the signal covariance matrix C ; c ™ b , m of the masked sky becomes 
non-diagonal in multipole space. Second, the noise level is not the same in different regions of the 
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sky, because the experiment looks at different regions for different amounts of time. This implies 
that also the noise covariance matrix Cf°^ l2m2 is non-diagonal in multipole space. The importance 
of the anisotropy of the noise was pointed out in jllj , as it caused an increase of the variance of the 
estimator £ at high Vs. 

For these reasons, it is worth trying to generalize the estimator £. We consider estimators for 
/nl (f° r a given shape of non-Gaussianities) which contain both trilinear and linear terms in the 
aim's. In the rotationally invariant case, a linear piece in the estimator can only be proportional to 
the monopole, which is unobservable. However, in the absence of rotational invariance, linear terms 
can be relevant. In this class, it is straightforward to prove that the unbiased estimator with the 
smallest variance is given by 

fiin(a) = ^'^2(( a hm 1 ai 2m2 ai am3 )iC limil4m4 C l2m2 ^ (23) 



N 

krrii 

—3 (ai imi ai 2m2 ai 3m3 )i (^ imi) j 2m2 C'j 3m3j j 4J7M oj 4m , 1 j • 
where N is a normalization factor: 

N = ^( a !imia| 2 m2 a l 3 m3)l ^dmi ,l 4 m 4 ^( 2 m 2 ,i 5 m 5 ^( 3 m 3 ,( 6 mfi ( a '4m 4 a l 6 m 5 a| 6 m 6 ) 1 • (24) 
knit 

The averages (. . .)i are taken with /nl = 1- As we stressed, the covariance matrix C = C cmh + C noisc 
contains the effects of the mask projection and the anisotropic noise and it is thus non-diagonal in 
multipole space. 

It turns out that the estimator £n n saturates the Cramer- Rao bound, i.e. it is not only optimal 
in its class, but it is also the minimum variance unbiased estimator among all the possible ones. To 
prove this, we expand the probability distribution in the limit of weak non-Gaussianity (which is 
surely a good approximation for the CMB) using the Edgeworth expansion |191 120| l2T]. 



P(a\fNh) = (1 - hL}Aai imi ai 2m2 ai 3m3 )i-z ^ ^ ====== . (25) 

~ da hmi da hrn2 da hm3 J J (2ir) N det(C) 



Now, following the same arguments of ^3] which considered the same problem in the rotationally 
invariant case, the estimator is optimal (in the sense that its variance saturates the Cramer-Rao 
bound), if and only if the following condition is satisfied: 

log -P(a | /NL) TT(t \(C („\ f \ 

77 = -F(/nl) (tiin(a) - /nl) , (26) 

d /nl 

where -F(/nl) is a generic function 3 of the parameter /nl- From the probability distribution (|2*3|l it 
is easy to check that the estimator £\ in is proportional to dlog P(a\f^i,) /d in the limit of small 
/nl- We conclude that £\\ n is an optimal estimator for a nearly Gaussian distribution. 

We now want to make some approximations to the optimal estimator £\\ n to make it numerically 
easier to evaluate. The full inversion of the covariance matrix is computationally rather cumbersome 
(although doable as the matrix is with good approximation block diagonal [H]). We therefore 

3 When cq. I|2()l) is satisfied, the function F turns out to be the Fisher information for the parameter /nl, 
F (/nl) = ((dlogP(a|/ NL )MNL) 2 ). 
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approximate C l a in the trilinear term of eq. (|23[) by masking out the sky before computing the 
ai^s and taking C as diagonal: 



where (a; imi a; 2m2 a; 3m3 )i is still given by the full sky expressions of the last Section. Once we have 
made this approximation for the trilinear term it is easy to prove that the choice for the linear term 
which minimizes the variance is 



Note that we must not approximate the covariance matrix C in the numerator as diagonal. This 
would leave only a term proportional to the monopole aoo (which is unobservable) , as in the rota- 
tionally invariant case. 

The normalization factor N is given by 



where / s k y is the fraction of the sky actually observed. As shown in [22], this correctly takes into 
account the reduction of data introduced by the mask for multipoles much higher than the inverse 
angular scale of the mask. The accuracy of this approximation has been checked on non-Gaussian 
simulations for the local shape in 

Let us try to understand qualitatively the effect of the linear correction. Take a large region of 
the sky that has been observed many times so that its noise level is low. This region will therefore 
have a small-scale power lower than average. Now, in a given realization depending on how the large 
scale modes look like, this long-wavelength modulation of the small-scale power spectrum may be 
"misinterpreted" by the trilinear estimator as a non-Gaussian signal. Indeed, for the local shape, 
most of the signal comes precisely from the correlation between long wavelength modes and the 
small scale power 0J. On the other hand, for equilateral shapes, as we noted, the signal is quite 
low on squeezed configurations so that we expect this effect to be small. The linear term measures 
the correlation between a given map and the anisotropies in the power spectrum, thus correcting for 
this spurious signal. Clearly the spurious correlation is zero on average, but the effect increases the 
variance of the estimator. 

We will apply the linear correction of the estimator only for the local shape, since, as we will 
verify later, it gives a very small effect for equilateral shapes. Following the same steps as in the last 
Section to factorize the trilinear term in the estimator, we get an explicit expression for the linear 
piece in the local case: 




(27) 




(28) 




(29) 




(30) 



where the two maps Sab(^, r) and SBB(n,r) are defined as 




(31) 
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S B B{n,r) = V ^4^^i,mi(^)%r^ y «2,m2(™)( a «imi««2m 2 ) • (32) 

km,i 

As discussed above, it is only the anisotropic part of the matrix (ai irni ai 2rri2 ) that gives a contribution 
to the linear part of the estimator. This matrix can be decomposed as 

( a h mi 0.l 2 m 2 ) = (lljmi^) + ( a "i mi a ?2°m2)' (33) 

where a^™ b are the a; m 's of a map generated with CMB signal only and then performing the mask 
projection, and af™ se are the ones associated with a map generated with noise only and then per- 
forming the mask projection. Both of these two matrices have an anisotropic component. For the 
map (cii^ a f^ 2 ) > this arises only because of the sky cut, while for the map (o"^^""™), it is gener- 
ated both by the sky cut and by the anisotropy of the noise power. We already discussed about the 
effect of the anisotropy of the noise in the previous paragraph; this effect turns out to be the most 
relevant. The sky cut gives a much smaller effect because, as we will explain more in detail in the 
next Section, it is mainly associated with the average value of the temperature outside of the sky 
cut, and this average value is subtracted out before the analysis. 



4 Analysis of WMAP 1-year data 

We keep our methodology quite similar to the one used by the WMAP collaboration for their 
analysis of the local shape in order to have a useful consistency check. We compare WMAP data 
with Gaussian Monte-Carlo realizations, which are used to estimate the variance of the estimator. 
We generate with HEALPix 4 a random CMB realization, with fixed cosmological parameters, at 
resolution nside=256 (786,432 pixels). The parameters are fixed to the WMAP best fit for a ACDM 
cosmology with power- law spectrum [21]: tt b h 2 = 0.024, VL m h 2 = 0.14, h = 0.72, r = 0.17, n s = 0.99. 
With these parameters, the present conformal time tq is 13.24 Gpc and the recombination time is 
tr = 0.27 Gpc. A given realization is smoothed with the WMAP window functions for the Ql, Q2, 
VI, V2, Wl, W2, W3 and W4 bands [23]. To each of these 8 maps we add an independent noise 
realization: for every pixel the noise is a Gaussian random variable with variance (Tq/N^, where 
Aobs is the number of observations of the pixel and do is a band dependent constant The maps 
are then combined to give a single map: we make a pixel by pixel average of the 8 maps weighted by 
the noise ofj/Aobs- We use this procedure because it is identical to that used in [2], thereby allowing 
direct comparisons between our results and those of the WMAP team. In future work, it can in 
principle be improved in two ways. First of all, for the window function to be strictly rather than 
approximately uniform across the sky, the weights of the eight input maps should be constant rather 
than variable from pixel to pixel. This is a very small effect in practice, since the eight A b s -maps 
are very similar, making the weights close to constant. Second, the sensitivity on small scales can be 
improved by using /-dependent weights [27]: for instance, at very high I, most of the weight should 
be given to the W bands, since they have the narrowest beam. This would also have a small effect 
for our particular application, since our estimator uses only the first few hundred multipoles. We 
explicitly checked that this would not reduce the variance of our estimator appreciably. 

4 See HEALPix website: http://www.eso.org/science/healpix/ 
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We apply the KpO mask to the final map, to cut out the Galactic plane and the known point 
sources ^Hl : this mask leaves the 76.8% of the pixels, / s k y = 0.768. The average temperature outside 
the mask is then subtracted. 

On the resulting map we calculate the estimators defined in the preceding Sections. For the local 
shape we performed the analysis both with and without the linear piece discussed in Section [3 The 
quadratic maps of equations (|31|) and (|32|) used for the linear correction are calculated by averaging 
many (~ 600) Monte-Carlo maps obtained with the same procedure above. We use HEALPix to 
generate and analyze maps at resolution level nside=256 (786,432 pixels). The integration over r is 
performed from tq — 0.025 tr up to tq — 2.5 tr with ~ 200 equidistant points, and then with another 
logarithmically spaced ~ 60 points up to the present epoch. Such a high resolution both in r spacing 
and in nside is necessary in order to reproduce the cancellation on squeezed triangles which occurs 
among the various terms in the equilateral shape l|14|l. The computation of each /nl on a 2.0 GHz 
Opteron processor with 2 GB of RAM takes about 60 minutes for the local shape, and 100 minutes 
for the equilateral shape. 

We then apply exactly the same procedure to WMAP data. These maps are analyzed after 
template foreground corrections are applied, in order to reduce foreground signal leaking outside the 
mask, as described in [T5| . 

A useful analytic bound on the variance of the estimators for /nl^ anc ^ /nl" 1 * s obtained from 
the variance of the full sky estimator (with homogeneous noise) (jSJ) with an / s k y -correction which 
takes into account the reduction in available information from not observing the whole sky. Taking 
into account the normalization factor, one readily obtains 

B 2 

-an 2 = /sky £ Q^fe ' (34) 

where -B; 1 z 2 z 3 must be evaluated for f^ff 1 or fm^ e( l ua l to 1- As we discussed, our approach is not 
strictly optimal as we are not inverting the full covariance matrix, so we expect o" an to be smaller 
than the actual standard deviation that we measure from our Monte-Carlos. 

In figure |1J we show the standard deviation of estimators for the local shape parameter 
as a function of the maximum multipole analyzed. We compare the results of our Monte-Carlo 
simulations (with and without the linear correction) with the analytic bound discussed above. The 
results without linear correction are compatible with the analysis in We see that the addition 
of the linear piece reduces the variance divergence at high Vs. The residual divergence is probably 
associated with the fact that we did not invert the full covariance matrix. The estimator with the 
smallest variance is the one with linear correction at Z max = 335, with a standard deviation of 37. 
The analytic approximation has an asymptotic value of 30 which should be considered the best 
possible limit with the present data. Our estimator thus extracts about (37/30) _2 fa 66% of the 
f^i^-uii brmation (inverse variance) from the WMAP data. The analysis of the data with l max = 335 
gives 47, so there is no evidence of deviation from pure Gaussianity and we can set a 2a limit of 

-27 < < 121 at 95% C.L. (35) 

In fig. EJ we show a map S^^n, r) for a radius around recombination calculated with noise only 
and for Z max = 370, as an illustrative example of the role of the linear piece. As we will clarify below, 
it is in fact the anisotropy of the noise that causes most of the contribution to the linear piece. The 
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Figure 4: Standard deviation for estimators of /^l as a function of the maximum I used in the 
analysis. Lower curve: lower bound deduced from the full sky variance. Lower data points: standard 
deviations for the trilinear + linear estimator. Upper data points: the same for the estimator without 
linear term, for which the divergence at high Z's caused by noise anisotropy had already been noticed 
in [Til- 



companion map SBB(n,r) is qualitatively very similar to the map Sab- It is the anisotropy of the 
noise that gives rise to a non-trivial (i.e. non-uniform) Sab map, as is clear from the comparison in 
the same figure with a plot of the number of observations per pixel. This can help us in understanding 
the contribution of the linear piece associated with the anisotropy of the noise to the estimator. Let 
us consider a particular Gaussian Monte Carlo realization with a long wavelength mode crossing 
one of the regions with a high number of observations. The trilinear piece of the estimator will then 
detect a spurious non-Gaussian signal associated with the correlation between this long wavelength 
mode and the small short scale power of the noise ( 5 ). The maps Sab and Sbb, because of their 
particular shape, will have a non-zero dot product with precisely the same long wavelength mode 
of the Monte Carlo map, and with an amplitude proportional to the anisotropy of the two point 
function of the noise, thus effectively subtracting the spurious signal detected by the trilinear piece 
and reducing the variance of the estimator. 

As the mask breaks rotational invariance, the CMB signal also gives non-uniform (S-maps. How- 
ever, the breaking of rotational invariance can be neglected far from the galaxy mask, so we expect 
that the contribution to the S'-maps from the CMB is to first approximation constant outside the 
mask (and zero inside). In this approximation, the linear contribution of the estimator would be 
sensitive only to the average value outside the mask. However, given that we are constraining the 

5 The effect of this spurious signal obviously averages to zero among many Monte Carlo realizations but it 
increases the variance of the estimator. 
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Figure 5: Top: SAB(n,r) map for r around recombination calculated with noise only. Bottom: 
number of observations as a function of the position in the sky for the Ql band (the plot is quite 
similar for the other bands). A lighter color indicates points which are observed many times. 

statistical properties of temperature fluctuations, the average temperature in the region outside the 
mask has to be subtracted before performing the analysis; therefore we expect the linear correction 
associated with the CMB signal to be rather small. These expectations are in fact verified by our 
simulations. We checked that the maps S coming from the CMB signal are to first approximation 
constant outside the mask, with additional features coming from the patches used to mask out bright 



15 




250 



200 



150 




100 



50 



250 



300 



350 



400 



450 



max 



Figure 6: Standard deviations for estimators of ' as a function of the maximum I used in 
the analysis. Lower curve: lower bound deduced from the full sky variance. Data points: standard 
deviations for the estimator without linear term. 

sources outside the Galactic plane. We then verified that if the average value of the temperature is 
not subtracted, the linear term coming from the CMB signal gives a very important reduction of the 
estimator variance, while its effect is almost negligible when the mean temperature is subtracted, as 
we do in the analysis. 

In figure EJ we study the standard deviation of the estimators for f^j^' ■ The estimator without 
linear corrections is quite close to the information-theoretical limit, so we did not add the linear 
corrections (which would require many averages analogous to the maps Sab and Sbb)- The reason 
why we have such a good behavior is that, as mentioned, here most of the signal comes from 
equilateral configurations. Thus the estimator is not very sensitive to the correlation between the 
short scale power and the (long wavelength) number of observations. Moreover the inversion of the 
covariance matrix is important only for the low multipoles, which give only a small fraction of the 
signal for f^ 11 '■ We find that the estimator with smallest standard deviation is at Z max = 405, with 
a standard deviation of 151. The analysis of the data with the same Z max gives —64. We deduce a 
2<7 limit 



The given limit is approximately 2 times stronger than the limit obtained in [I], indirectly obtained 
starting from the limit on /m 3,1 . The limit appears weaker than for the local case because, for the 
same /nl> the local distribution has larger signal to noise ratio than an equilateral one, as evident 
from fig. This is not a physical difference but merely a consequence of our defining /nl at the 
equilateral configuration. 

To check that the two limits correspond approximately to the same "level of non-Gaussianity" 



366 < /; 



•equil. 



< 238 at 95% C.L. 



(36) 
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we can define a quantity which is sensitive to the 3-point function integrated over all possible shapes 
of the triangle in momentum space. In this way it will be independent of which point we choose for 
normalization. We define this quantity, which we will call NG, directly in 3d as an integral of the 
square of the functions in fig. ^ 

NC _ ( f x 2 dx 2 x 3 dx 3 F(l,x 2 ,x 3 ) 2 
[J A (2-) 2 A|x 2 ^3 3 

The integration is restricted to the same triangular region as in fig.s Handel 1 — x^ < x^ < x^ < 1. 
The measure of integration comes from the change of variables from the 3d wavevectors to the ratios 
x 2 = k 2 jk\ and x 3 = k 3 jk\. NG is parametrically of order A^ 2 ■ /nl (where A<j> ~ 1.9 • 10~ 8 
|241 I28| ) , which is the correct order of magnitude of the non-Gaussianity, analogous to the skewness 
for a single random variable. To better understand the meaning of the defined quantity, note that 
if we integrate inside the parentheses in eq. (|37|) over the remaining wavevector we get the total 
signal to noise ratio (or better the non-Gaussian to Gaussian ratio). This further integration would 
approximately multiply NG by the square root of the number of data. This means that to detect 
a given value of NG we need a number of data of order NG -2 , in order to have a signal to noise 
ratio of order 1. This is clearly a good way to quantify the deviation of the statistics from pure 
Gaussianity. The 2a windows for /nl can be converted to constraints on NG (we define NG to have 
the same sign of /nl) 6 

-0.006 < NG local < 0.025 at 95% C.L. (38) 




-0.016 < NG cquiL < 0.010 at 95% C.L. (39) 

Contrary to what one might think from a naive look at the limits on /nl> the maximum tolerated 
amount of non-Gaussian signal in a map is thus very similar for both shapes. 

So far in this paper, we have neglected any possible dependence on the cosmological parameters. 
A proper analysis would marginalize over our uncertainties in parameters and this would increase 
the allowed range of /nl- In order to estimate what the effect of these uncertainties is, let us imagine 
that the real cosmological parameters are not exactly equal to the best fit ones. The cosmological 
parameters are determined mainly from the 2-point function of the CMB, C; cmb , therefore the largest 
error bars are associated to those combinations of parameters which leave Cf mh unchanged. In the 
limit in which the Cf mb, s remain the same, also the variance of our estimator (which is computed 
with the best fit cosmological parameters for a ACDM cosmology with power-law spectrum |24j ) 
remains the same: in the weak non-Gaussian limit the variance just depends on the 2-point function. 
However, the combination of parameters which leave unchanged the Cf mb, s does not necessarily 
leave unchanged the bispectrum, therefore uncertainties in the determination of the cosmological 
parameters will have an influence on the expectation value of the estimator. The expectation value 
of the estimator would be: _ 

< £ >4 £ • (40) 

ll<l2<h 11 12 13 

6 Note that for the local model the integration in eq. 1)37(1 diverges like the log of the ratio of the minimum 
to the maximum scales in the integral. For the quoted numbers we chose a ratio of 1000. 
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where N, -B; 1 ; 2 z 3 and Ci are the same as in our analysis, computed with the best fit cosmological 
parameters, while -Bz 1 ; 2 ; 3 is the true bispectrum. Thus when varying parameters, the normalization 
N needs to be changed to make the estimator unbiased. 

The most relevant uncertainty is the reionization optical depth r, which is correlated with the 
uncertainty in the amplitude of the power spectrum A<j>. With the purpose of having a rough 
estimate, we can approximate the effect of reionization by a multiplicative factor e~ T in front of the 
transfer function A^(fc) for I corresponding to scales smaller than the horizon at reionization. In 
order to keep the Cf mb 's unchanged, at least at high i's, we then multiply A$ by e 2r . The bispectrum 
is proportional to A|, Af(k) 3 and thus changes even if the Cj's do not. For the equilateral shape 
most of the signal comes from equilateral configurations with all the 3 modes inside the horizon at 
reionization; in this case (£) scales roughly as e T . For the local shape the signal comes from squeezed 
configurations with one mode of much smaller wavelength than the others. Taking only 2 modes 
inside the horizon at reionization we get (£) oc e 2r . 

If we consider the la error, r = 0.166lQgy ( j i [21], we find that the uncertainty in the reionization 
depth should correspond to an error on /nl of order 8% for the equilateral shape, and 15% for the 
local shape. These numbers translate directly into uncertainties in the allowed ranges we quote. 
In the future, if the error induced by the uncertainty in the cosmological parameters will become 
comparable to the variance of the estimator, a more detailed analysis will be required; however, at 
the moment an analysis with fixed cosmological parameters is certainly good enough. 

5 Conclusions 

We have developed a method to constrain the level of non-Gaussianity when the induced 3-point 
function is of the "equilateral" type. We showed that the induced shape of the 3-point function 
can be very well approximated by a factorizable form making the analysis practical. Applying our 
technique to the WMAP first year data we obtained 

-366 < /° q L uil - < 238 at 95% C.L. (41) 

The natural expectation for this amplitude for ghost or DBI inflation is /nl 1 ' 1 ~ 100) below the 
current constraints but at a level that should be attainable in the future. The limit an experiment 
can set on /nl just scales as the inverse of the maximum I the experiment can detect. As a result, 
in the case of WMAP, increased observing time will approximately decrease the error bars by 30% 
and 60% for 4 years and 8 years of observation. The increased angular resolution and smaller noise 
of Planck pushes the point where noise dominates over signal to I ~ 1500. This should result in 
a factor of 4 improvement on the present constraints. In addition polarization measurements by 
Planck can further reduce the range by an additional factor of 1.6 29 . 

We also constrained the presence of a 3-point function of the "local" type, predicted for example 
by the curvaton and variable decay width inflation models, obtaining 

-27 < ,/^ al < 121 at 95% C.L. (42) 

We defined a quantity NG, which quantifies for any shape the level of non-Gaussianity of the 
data, analogously to the skewness for a single random variable. The limits on NG are very similar 
for the two shapes, approximately |NG| < 0.02 at 95% C.L. . 
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We showed that unless one has a full sky map with uniform noise, the estimator must contain 
a piece that is linear on the data in order to extract all the relevant information from the data and 
saturate the Cramer-Rao bound for the 3-point function measurement uncertainty. This correction 
is particularly important for the local shape and accounts for the improvement of our limits over that 
from the WMAP team. Moreover this correction goes a long way towards reducing the divergence 
in the variance of the estimator as Z max is increased. 
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